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ABSTRACT 

This paper presents a model which can be used to predict the response 
of wind turbines to atmospheric turbulence. The model was developed 
using linearized aerodynamics for a three-bladed rotor and accounts 
for three turbulent velocity components as well as velocity gradients 
across the rotor disk. Typical response power spectral densities are 
shown. The system response depends critically on three wind and turbu- 
lence parameters, and models are presented to predict desired response 
statistics. An equation error method, which can be used to estimate the 
required parameters from field data, is also presented. 


WIND TURBINE SYSTEM MODEL 

Before embarking on a discussion of the deti.5 ed characteristics of 
atmospheric turbulence pareuneters, it is necessary to present the model- 
ing framework in which the parameters will ba used to predict system 
responses. The primary purpose of the model is to provide a tool by 
which designers can estimate the effects of fluctuating turbulence 
inputs on the wind turbine, structural and power system responses. 

For an n degree of freedom system, the basic principles of Newtonian 

mechanics tl] give equations of motion of the form 

tM]{z} + tC Hz} + [K ]{z} = {f } (1) 

s s a 

where {z} = the nxl vector of generalized displacement coordinates* 

[M] = the nxn inertia matrix. 

[C ] = the nxn gyroscopic and structural and power train damp- 
^ ing matrix. 
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[K ] « the nxn structural and power train stiffness matrix, 
{f®} » the nxl vector of aerodynamic forces and moments gener- 
ated by the turbine rotor . 


The aerodynamic forcing term of Eq. (1) depends upon the motion of the 
turbine rotor with respect to the ground as well as the motion of the 
air. If the aerodynamic forces and moments are linearized about a 
steady operating condition, the following equation results 


{f } = {f } + [F]{u} - tc Hz} - [K ]{z} 
an a a 


( 2 ) 


where 


{f } = the nxl vector of steady, nominal aerodynamic forces 
and moments. 

{u} “ the mxl vector of fluctuating turbulence Inputs. 

[F] = the nxm matrix of aerodynamic influence coefficients. 
[C ] = the nxn aerodynamic damping matrix. 

[K*] « the nxn aerodynamic stiffness matrix. 


In this particular model, the turbulence input vector {u} consists of 
three velocity components which are uniform over the turbine rotor disk 
and sik additional gradient terms which account for variations in tur- 
bulent velocity over the rotor disk. Tedale 1 gives a verbal description 
of the nine turbulence input terms appropriate for a rigid, three-bladed 
wind turbine rotor. 


TABLE 1. DESCRIPTION OF TURBULENCE INPUT TERMS 


Component 


Description 


V 

X 

V 

y 


V 

z 

V 

y»x 

V 



uniform lateral or side component (in rotor plane) 

uniform longitudinal component along steady wind 
direction 

uniform vertical component (in plane) 
lateral gradient of longitudinal velocity 
vertical gradient of longitudinal velocity 
swirl about steady wind axis (in plane) 

shear strain rates (in plane) expressed in a refer- 
ence frame rotating at three times the rotor rate 

in-plane dilation 
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Atinuminq that tho atmoaphorlc' tuibuXunca li, aduquatolY a.Miurihod by th.’ 
homoqunoouo, iaot-iropic Von Karman modol i:’l, 

can Lh) aptn-oxltnatod by tho foUowinii not ol ntochantic airfounUal 
uquat.ionii IXl 


lul •• lA^llu.1 + lb Hwl 
w ^ 


(d) 


whoro 


Iw} 


i». 


an mxl. vootor ot whito luriao oxcitutloiiH with tlat 
ixiwor apoctral donsity, B " 

the mxm dynamica matrix for tho turbulonoi' input a. 
tho mxm distribution matrix for tho whlto noise oxci- 
tations. 

Tho matrices IA„1 and IB„] are diagonal, except for two off diagonal 
tvZ il l/Vl, which account for the throe rotations per rotor rovol.- 
tion effect in tho Cj. and Yj. terms caused by tho throe blades moving 
through the in-plano turbulence gradients. 

The motion Kgs. (1), the aerodynamic force Kgs. (2), and tho wind turbu- 
lence inputs Kg«. (i) combined into a sot ot system oguatxons 

of tho form 


whort 


lx) - IMixl + iwHwl 


{y} --- icUx) + IYj^> 


(4) 


Ixl - 
Iwl 

ly] = 

lv„> 

lAl 


Ihl 


tho Nxl t?yi:itom stato vootor (N ~ Jntm) 

^tho n\xl white noise turbulence excitation vector, 
tho Qxl vector of system ros)xinso variables, 
the <xl vector of steady nominal system responses. 


0 

1-M*'(K„+K.) -m"(C to ) 

O'’ lO' 


0 

m*'k 

A 

W 


tho NXN e^y^Hom 
matrix 


w J 


tho Nxm whito noii^o oxoitiU ion diiU r ibut ion 
matrix* 


ICl 


tho V'xN roTiponfvO dialr ibut io!\ matrix. 


Note tlut 0 !r. and are deviat ions from the steady, generalised dis- 
placement and velocity comixments. The ouli'uts lyl a.ui the correspo.u - 
ing matrix ICl depend upon the particular s.'t ol displacements, vcKh i 
ties or load response vai ial'les of interest to the desianei. 

Till' i’.ystem I'lpiat ions, of mot ion aiven by l.ii. t'i' aii di i i\< d .i...>umini 
rigid, t hree-bl .ided turbine voior. It is. pos.sible to de. ive s.ys.t em 
eguations tor two-bl.ided rotors similar to t hes.e eguat ions , except 
th.it several of the terms in the |Al matrix will li.ive periodic terms 
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instead of being constant as in Eq. (4) . 


At this point, wo will describe briefly how the wind paroroetera ontor 
the various coefficients of the overall system model. First, the 
steady wind speed, Vw/ affects the nominal aerodynamic forces and the 
linearized aerodynamic coefficients in the matrices iCg] , iKaJ and (F] . 
Second, both the steady wind speed, and the turbulence integral 
scale, L, affect the matrices lA^,] and . Finally, the turbulence 
component variance, c^, as well as Vw and L, affect the power spectral 
density, S^,, for each of the white noise excitation components. Thus, 
three atmospheric turbulence parameters, v„, c, and L, must be known 
in order to utilize the model given by Eq. (4) . 


Once the appropriate turbulence parameters are specified, the response, 
power spectral densities can be computed using the model given by 
Eq. (4) . Since the white noise inputs are uncorrelated, the following 
equation results 


{S^(w)} = [T(u)Hs^} 


(5) 


where 


{S (w) } 
ITloj)] 


the )lxl vector of response power spectral densities 
the mxl vector of white noise excitation power 
spectral densities. 

the £.xm matrix of squared, complex magnitudes of 
the system frequency response matrix elements, 
the radian frequency. 


If 


is one element of [T((o)], then 




( 6 ) 


where H,. (ito) « the corresponding element of the complex frequency 
^ resp onse matrix, 

i = 


Assxaning the eigenvalues of the system dynamics matrix lA] are distinct 
the complex frequency response matrix is given by 


tH(iu)] = tC] IM] [iwtil - tA]f'[Mr'tB] 


(7) 


where 


IM] = complex modal matrix consisting of columns of eigen- 
vectors of lA] . 

[A] = diagonal complex matrix of eigenvalues of lA] . 

[I] identity matrix. 
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TYPICAL WIND TURBINE RESPONSE CHARACTERISTICS 

A siroplifiod fivQ-dQgrae-of-froQdom model for a throe-bladod, horizontal 
axis wind turhino was dovoloped by Thresher > ot al, 14]. The five 
gonoralized displacement degrees of freedom are given by 


{z}’' = 


(8) 


where U 
V 

<l> 

X 


■=• lateral displacement of the nacelle in x direction. 

= fore-aft displacement of the nacelle in y direction. 
= yaw angle. 

= pitch angle. 

=■ rotor angular displacement about spin axis. 


Figure 1 shows the coordinate system used for this model. The configu- 
ration shown in the figure is appropriate for a down-wind rotor design. 


Thresher and Holley [5] utilized the model for two typical wind turbines 
of widely differing size. The first, designated the Mod-M, is an 8 kW 
free yaw system with a down-wind rotor. The second, the Mod-G, is a 
large 2.5 MW machine with a fixed yaw, up-wind rotor. The system 
characteristics for these two machines are given in Tables 2 and 3. 


TABLE 2. MOD-M CHARACTERISTICS 


Rotor Characteristics ; 

Rotor Radius 
Hub Height 

Blade Chord (constant) 

Coning Angle 

Blade Twist 

Pitch Setting (to ZLL) 

Steady Operating Conditions ; 

Rotor Speed, D 
Wind Speed, 

Approximate Output 

Aerodynamic Properties > 

Lift Curve Slope 
Drag Coefficient Cp^ 

Turbulence Parameters ; 

Standard Deviation, o 
Integral Length Scale, L 

System Frequencies (Tower Motion) ; 

1st Bending (fore-aft) 

2nd Bonding (fore-aft) 

1st Bending (side-to-side) 

Ist Torsion 


5.081 

m 

(16.67 ft) 

16.8 

m 

(55 ft) 

.457 

m 

( 1.5 ft) 

.061 

rad 

( 3.5‘») 

0 

rad 

( 0.0») 

.052 

rad 

( 3.0®) 


7.681 rad/s 

(73.35 

RPM) 

7.434 m/s 

(16.63 

MPH) 


6 kW 


5.7 

.02 


.619 

m/s 

( 2.03 ft/s) 

91.44 

m 

(300 ft ) 

15.1 

rad/s 

(2.0 fi) 

53.1 

rad/s 

(7.0 il) 

15.9 

rad/s 

(2.1 fi) 

0.0 

rad/s 

(Free Yaw) 
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TABLE 3. MOD-G CHARACTERISTICS 


Rotor CharactorifiticB> 


Rotor Radius 

45,7 m 

(150 

ft) 

Hub Hoight 

61.0 m 

(200 

ft) 

Blade Chord (linear taper) 

2.36 m 

( 7.74 

ft 


to . 96 m 

to 3.15 

ft) 

Coning Anglo 

.070 rad 

( 4.0") 


Blade Twist (linoar) 

.140 rad 

( a.o“) 


Pitch Setting at Tip (to ZLL) 

.108 rad 

( -6.2») 


Steady Operating Conditions* 

Rotor Speed, fi 

1.833 rad/s 

(17.5 

RPM) 

Wind Speed, 
Approximate Output 

8.940 m/s 
1.1 MW 

( 20.0 

MPH) 

Aerodynamic Properties t 

Lift Curve Slope 
Drag Coefficient, Cp^ 

5.73 

.008 



Turbulence Parameters: 

Standard Deviation, a 

.744 m/s 

• 

CM 

ft/s) 

Integral Length Scale, L 

152.4 m 

(500 

ft ) 

System Frequencies (Tower Motion) : 

1st Bending (fore-aft) 

2.75 rad/s 

( 1.5 

«) 

2nd Bending (fore-aft) 

12.8 rad/s 

( 7.0 

«) 

1st Bending (side-to-side) 

2 . 9 rad/s 

( 1.6 

ft) 

1st Torsion 

9.5 rad/s 

( 5.2 

ft) 


Two aerodynamic wake models were used for each system to compute the 
coefficients in the aerodynamic system matrices [C^l # and [P] . 

In the first, the steady conditions are used with standard momentum 
theory to compute the steady distribution of induced velocity across 
the rotor disk. This induced velocity is then assumed constant for 
the given conditions. This model is called the "Frozen Wake." In the 
second model, the induced velocity which results from a slowly varying 
velocity field is computed using a quasi-steady momentum balance. In 
this model, the turbine thrust is always in equilibrium with the driving 
turbulent velocity, and is called the "Equilibrium Wake." Aerodynamic 
stall is not modeled in either caL.e. 

Figures 2 and 3 show the power spectral densities of the thrust load 
and the yaw angle for the Mod-M machine. In the low frequency portion 
of Figure 2, the thrust load response closely follows the power spec- 
trum of the Vy turbulence input. At higher frequencies the resonance 
effects of the tower bending modes are observed. In Figure 3, the yaw 
response is dominated by the turbulence input. This turbulence 

input term can be interpreted as the rate of change of the direction 
in the horizontal turbulent velocity component. A smaller additional 
effect is duo to the uniform side velocity, V^, turbulence term. 
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Figure 3. Yaw Response, (j), for Mod-M. 
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PigmroR 4 an4 5 show airoilai: rooulta for tho Mod^G maohino oxoopt that 
thn towor torsion load is shown instead of tho yaw angle for this fixed 
yaw maohino. The Mod-G maohino shows a greater sensitivity to tho 3U 
off eats of tho e and y turbulonao terms # 


MBTHOPOIiOGV FOR COMPUTATION OP RESPONSE STATISTICS 


This sootion gives a brief disounoion of tho toohniquon by which tho 
model given by Eq. (4) oan bo used to oomputo dooirod rooponoo stotio” 
tics. Assuming that tho fluctuating components of tho atmospheric 
turbulence are adequately described by Gaussian statistics [0^7] i the 
model will give the conditional probability density function of tho 
response given tho steady wind speed Vw, and tho turbulence parameters 
o and L. Thus^ considering only a single, scalar response variable 


p(y|v^,o,L) 


1 -l/2(y-M /o 

— O ^ ^ 


(9) 


where 


M ® y (V ) “ the steady response for given V 
y n w w 

Cy « Oy(V^ a,L) “ the rms response for given V^# o, and L. 


This function can be recognized as the standard Gaussian density func- 
tion. The conditional mean, Vy, is a nonlinear function of V^, and 
the conditional rms response, Cy, depends nonlinear ly on Vw and L and 
is proportional to a. The rms response, o^, can be computed from the 
response power spectral density by the relation 


o 


2 

y 


1 

TT 


w 

/ S (t»j)d(0 
o ^ 


( 10 ) 


The response variance o 
latlon [8] 


2 

y 


can also be calculated directly using the re- 


cTy = tCHMHP] 


( 11 ) 


where 


to “ thu row matrix relating the response to the system 
state vector. 

(M] “ modal matrix with column eigenvectors. 

* « complex conjugate of the matrix. 


The Hermitian matrix, [P] , satijf ies the linear relation 


where 


(AHPl + tPHA]* + ImF’ (“) - 0 

w 

(A) = the diagonal matrix of complex eigenvalues. 

(B) » tho white noise input distribution matrix. 


( 12 ) 
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Notu (hat tho matrloon 1MI» |M1 , mul [A] tlapontl nonllnaaily on th*> 

tKirami'l-orn V aiul li* 

* w 

Now, mM>i>otu« H lit aotiin'il to c-omputo Mu' probability that, y ttxooodn a 
o.'rl-aln orlt loa.! valno y^,. Ttio ooncUtional proUabiUty is qlvon by 


rrty ^ y Iv ,o,l.l " / pCyjv ,o,b)dy 
o w w 

Hubiitltul Inq Hij. (") Into Kq. (13) yloUtn 


whoro 


Priy 


Y,JV, 


,a,h 


w 


y 

orf (— 


c 



ort’(») 


1 



(•) 

/ 



o 


the error function. 


(. 13 ) 


( 14 ) 


The total probability ia thua qivon by 


tv ix' w y 

Pr { Y 1 y . ) = / / / (7 “ f (-^—-^) ) V (V^, 0 , b) dV^dodL ( lb ) 

000“ Y 

whori' p(V ,o,L) - the ioint probability density function of the 

ivsitive wind and turbulence parameters. 

Por oomput at lonal puriKtses, the inteqrals can be approximated by dis- 
crete aunmui lions, so that 


Prly y . } 




(” - orf(— 

y 


) )l’^Vw^fnj^,L^j) 


( 16 ) 


where tlu' subscripts denote discrete Viilues of the t’^irametors associated 
with "countinq bins.” The probability required is the joint probability 

that V is in bin j, 0 is in bin k, and L is in bin >! . 
w 

Ui\l ort.unaioly , complolo data fur dot urmininsi tlio Joint denaity funct ion 
ti>r I ho wind and turbuUmco paranu'lt'ra in qonorally laokinq. Howovor, 
Movt'val aimplifyimi anMumptiona make an iU>proximato modol iK^saiblo* 

in iin iitmoaplii'r io lunuidary layor with nouiral buoyant atability llu> 

Uniiii iihmlo proi ilo haa boon found to adoquali'ly modol the variation 

of V wtlli lu'iithi hM . Thin modol i it of Iho form 
w 


wlu'rc' u^ triotik'n 

;• hoiqht abovo t ho around, 

;• iuMuit\al hoiqhl vfhcrv V 0 (i>tl».'n /,<'ro) • 


oniGiriAL H 

» terrain roughness length. POOR QU/\Un 3 

Pjjost; at al. [10] recommend the Welbull probability distribution for 
the steady wind speed at the reference height of 10 m. Thus solving 
Eq. (17) for u^ when z = = 10 m yields 

z-z +z 


w 


- . no. 
in( — ::: ) 


z -z +z 
. . r n o. 
in ( ) 


(18) 


where 


V at the reference height, 
w 

reference height. 


Since V„ and V- are linearly related, also satisfies the Weibull 
distribution which can be differentiated to give the density function 


of the form 


k ,\k-i -<w 

p<V ■ - ^ 

o o 


(19) 


where 


» a site parameter (“ 2) . 


w 


w 


r(l + ^) . ^ 

annual mean wind speed at the desired height. 


r(») = gamma function. 

The annual mean wind speed at the desired height can be found from the 
value at the reference height by the use of Eq. (18) . 

The rms, turbulent component velocity, o, is found to be highly corre- 
lated with the steady wind speed. Panofsky, et al. [11] give the re- 
lation 


o = 2.3 u^ 

so that when Eq. (17) is used for u^. 


( 20 ) 


a 


0.92 


z-z +z 
- - r o, 

^■»^( — r > 


V 

w 


(21) 


The turbulence integral scale, L, is much less understood. Most evi- 
dence indicates that it is independent from the steady wind speed, V^, 
and the variance, o^. Several authors [12,13,14] recoiranend different 
power laws for the variation of integral scale with height. However, 
these relations are inconsistent and the experimental data exhibit 
wide scatter. It is highly recommended that an experimental program be 
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undertaken to determine an appropriate height scaling law and to account 
statistically for the variation observed at a given height. In the 
interim, we will assxame the integral scale is deterministic and satisfies 
the height relation 


L » L 


r 



where 


=> a site parameter (.“ 65 m) . 
« reference height <=• 10 m. 


( 22 ) 


Using these simplifying approximations for the parameter models, the 
statistical procedure given by Eq. (15) reduces to 


Pr{y ly } » / (| - erf (-2-^) ) p (V^) dV^ (23) 

° o y 

where given by Eq. (19). 

The quantities y and o will be complicated functions of V„ given by 
the model of the^ turbine response, with Eqs. (21) and (22) used for the 
parameters a and L. Obviously, numerical procedures would be used to 
perform this computation. 


ESTIMATION OF MODEL PARAMETERS PROM FIELD DATA 

Since the steady wind and turbulence parameters, V„, a, and L, criti- 
cally affect the statistics of the response, it is highly desirable to 
have a reliable method for extracting the parameters from real field 
data. One such method is the equation error method [15]. Basically, 
the method determines a set of parameter values which minimize the 
difference between the data and predicted values based on the model 
equations. The resulting parameters will then serve to characterize the 
turbulence sample observed. A whole collection of such parameter values 
will then give the required statistical information discussed in the 
previous section. 

Before proceeding to give the detailed procedure for estimating the mean 
wind and turbulence parameters, a brief description of the equation 
error method will be given. Suppose we have an accurate, noise-free 
measurement of a random process, u, modeled by the stochastic differen- 
tial equation. 


u = au + bw 


(24) 


where w = white noise with flat PSD = S^. 
a,b = model parameters. 
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The measurements will be a set of N values^ u(i) taken at discrete 
times with a constant time interval, x, between measurements. The 
continuous time model can be converted to the discrete time form 


u(i+l) = e®”^ u(i) + ?(i) 


( 25 ) 


where 5(i) = a random sequence of uncorrelated values. 

2 

The variance of 5(i) is found by matching the stationary variance of 
u(i) andu(t). Thus, from Eq. (25) 

Elu (i+l)3 = e'^®^ Elu (i)] + Eir(i)] (26) 


which when solved yields 


= El^^(i)] 


,, 2ax. 2 
(1 - e )o^ 


(27) 


Prom Eq. (24) (assuming a < 0) , 


2 2 

2a<j + b S “0 
u w 


Using Eq. (28) in Eq. (27) yields 

2 ,, 2at. , b^ _ . 

= (1 - e ) (- - S^) 


(28) 


(29) 


Now, since u(i+l) and u(i) are linearly related and the noise term is 
sequentially uncorrelated, standard regression methods [16] can be 
used to estimate e®T and cr| from the data sequence. Thus, we choose 
the parameter, a, to minimize the estimated variance 


"2 


1 

N-1 


The product. 



2 

b S 

w 



• Z (u(i+l) - e^'^u(i)) 
i=l 

is determined from Eq. 
■'2 

-2a a 



2a X 

- e 


2 

(29) 


(30) 


(31) 


It is impossible to estimate b and separately. 

With the mathematical preliminaries out of the way, let us return to 
the turbulence parameter estimation problem. Suppose wo have two 
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propeller typo anemometers set up to measure orthogonal horizontal 
components of the wind. Let Vj^(i) and V 2 (i) be sequences of measure- 
ments taken from the anemometers. The first step in the procedure is to 
find the steady wind speed and direction. Thus, determine 


<Vi> 



N 

I V (i) 
i=l 


<V2> 


1 

N 


N 

I 

i-1 


v^Ci) 


(32) 


Now, 


w 


} 2 
/<v, > + 


ij) = tan 


'r • "^2" 
-1 "^2> 


<vi> 


(33) 


The lateral and longitudinal turbulence components are thus determined 
from 


V (i) => v^(i) cos(|) - V, (i) sinij) 

X 2 1 

V (i) = V (i) cos(j) + v_(i) sin«|) - V 

y X ^ w 


(34) 


The next step is to determine the parameter, L, using the equation error 
regression procedure. According to the model developed by Holley [17] , 
the lateral and longitudinal components of the turbulence satisfy the 
stochastic differential equations 




V 

y 



/2 


w 


w. 


(35) 


where Wj^ and W 2 are independent white noise processes with equal power 

spectral densities, S = o^L/V^. 

w w 


Applying the equation error regression technique of Eq. (30) and normal- 
izing each of the equation errors by the variance gives the variance 
estimate 


"2 

a 


1 

2 


"2 "2 

^1 ^2 

-4V t/L -2V t/L^ 

w , w 


(36) 
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where 


-^2 1 

o: = ~r 2 (V (i+1) - e 

^ i»l 

.2 1 N-1 


"2VX/L 2 


-V t/L 2 

” V (1))^ 

y 


The value of L is chosen to minimize a and a is the resulting a after 
the minimization. 


The parameter values determined by this method will characterize ^e 
particular turbulence sample observed during a given sampling period. 

It is expected that the values will be different for different days and 
times at which the data is taken. This collection of parameter values 
can then be used to estimate the statistics discussed in the previous 
section. 


CONCLUPIONS 

The paper has presented a modeling technique which can be used to esti- 
mate wind turbine response statistics due to atmospheric turbulence. 

Up to this point all of the modeling results have been theoretical. 
Before these techniques can be put to use by designers, it is required 
that they be verified using atmospheric and wind turbine field data. 
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